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ABSTRACT 

The  numerical  synthesis  of  regional  seismograms  has  become  an  integral  part  of  Lawrence 
Livermore  National  Laboratory’s  (LLNL)  seismic  discrimination  program.  In  this  paper,  we 
summarize  our  fundamental  approaches  to  numerical  modeling.  Our  capabilities  currently 
include  reflectivity,  normal  mode,  boundary  integral,  and  finite-difference  modeling,  along  with 
hybrid  approaches  which  utihze  two  or  more  of  these  techniques  together.  We  apply  these 
capabilities  to  the  discriminant  variability  along  three  different  arrays  deployed  during  the  Non- 
Proliferation  Experiment  (NPE).  Phase  amplitudes  have  been  calculated  for  the  approximately 
three  hundred  regional  stations  which  recorded  the  NPE.  The  majority  of  these  recordings  were 
to  the  west  of  the  NPE,  along  one  profile  of  the  Southern  Sierra  Continental  Dynamics  (SSCD) 
refraction  profile.  Based  on  the  three  profiles  which  made  up  the  SSCD  refraction  experiment, 
traveltime  tomography  was  utilized  to  develop  a  well  constrained  3D  velocity  model  across  a 
profile  which  extends  from  the  Basin  and  Range,  through  the  southern  Sierra  Nevada  range,  the 
Great  Valley,  to  the  San  Andreas  fault  zone  in  the  coastal  ranges.  The  western  array  of  the  NPE 
regional  deployment  consisted  of  285  stations  coincident  with  this  SSCD  profile.  This  resulted 
in  phase  and  discriminant  coverage  along  one  of  the  most  well  constrained  velocity  profiles  in 
the  western  United  states.  Our  analysis  shows  that  although  there  is  amplification  of  Pg  and  Lg, 
Pn  has  some  of  the  most  dramatic  variations  in  amplitude.  In  California,  these  amplifications 
coincide  with  the  western  flank  of  the  Sierra  Nevada,  the  eastern  quarter  of  the  Great  Valley, 
and  the  coastal  ranges.  In  Nevada,  a  dramatic  amplification  occurred  in  the  Broken  Hills 
volcanic  region.  To  the  east,  along  the  Arizona-Utah  border,  amplification  occurred  along  the 
transition  between  the  Basin  and  Range  and  the  Colorado  Plateau.  The  amplification  of  Pn  in 
these  regions  resulted  in  a  considerable  deviation  in  6-8  Hz  Pn/Lg  ratios.  However,  we  find  that 
the  Pn/Lg  slope  was  considerably  more  stable,  showing  less  variation  as  a  function  of  distance 
from  the  NPE.  Numerical  modeling  with  the  reflectivity  shows  that  when  material  contrasts  are 
large  enough,  a  local  plane  layered  structure  can  result  in  significant  Pn  amplification.  Boundary 
integral  modeling  shows  that  site  resonances  due  to  local  2D  structure  and  topography  can 
account  for  the  dramatic  amplification  and  time  duration  of  Pn  observed  in  Broken  Hills.  Finite- 
difference  modeling  shows  that  crustal  thickening  along  the  Colorado  Plateau  can  result  in 
focusing  which  may  explain  much  of  the  amplification  in  the  0.5  to  10  Hz  frequency  band  which 
appears  at  distances  of  200  and  280  km. 
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OBJECTIVES 

A  primary  goal  of  our  research  at  LLNL  has  been  to  develop  an  assemblage  of  the  most  promising 
numerical  techniques  for  simulating  regional  wave  propagation  through  complex  media.  These  numerical 
capabilities  will  be  an  integral  part  of  our  approach  to  regional  characterization  in  both  the  Middle  East 
and  North  Africa.  In  this  paper,  we  focus  on  demonstrating  and  validating  our  current  modeling 
capabilities,  which  include  reflectivity,  boundary  integral,  finite-difference,  along  with  hybrid  forms  of 
these  approaches.  We  use  the  boundaiy  integral  and  finite-difference  techniques  to  provide  better 
understanding  of  regional  discriminant  variability  along  three  arrays  which  were  deployed  to  record 
regional  signals  from  the  NPE.  These  arrays  extend  west,  northwest,  and  east  from  NTS,  with  the 
western  line  coinciding  with  one  of  three  Southern  Sierra  Continental  Dynamics  (SSCD)  refraction 
experiment  profiles  (Ruppert  et  al.,  in  press;  Fliedner  et  al.,  in  press).  This  gives  detailed  coverage  of 
regional  phases  along  one  of  the  most  well  constrained  crustal  profiles  in  the  western  United  States. 
Phase  amplitudes  and  the  resulting  discriminant  variability  are  presented.  This  “ground  truth”  dataset  will 
then  act  as  the  basis  for  the  future  validation  of  our  numerical  codes  concurrent  with  our  regionalization 
of  the  Middle  East  and  North  Africa. 

RESEARCH  ACCOMPLISHED 

Current  Capabilities 

Synthetic  modeling  of  seismic  wave  propagation  in  elastic  media  is  an  integral  part  of  LLNL's  seismic 
discrimination  program.  The  purpose  of  our  research  is  to  apply  a  collection  of  validated  wave 
propagation  codes  to  the  Middle  East  and  Northern  Africa  for  discriminant  phase  analysis.  This  includes 
discriminant  development,  the  interpretation  of  discriminant  results,  and  regional  characterization.  Our 
current  numerical  capabilities  for  seismic  wave  propagation  include  a  generalized  ID  reflectivity 
algorithm  (Randall,  1994)  and  a  2D  boundary  integral  algorithm,  which  generalizes  the  reflectivity 
approach  to  irregular  boundaries  (Schultz  and  Toksoz,  1994;  Bouchon  et  al.  1989).  Our  current  finite- 
difference  capabilities  consist  of  a  highly  optimized  2D/3D  wave  propagation  code,  ELAS3D  (Levander, 
1988;  Larsen  and  Harris,  1993;  Larsen,  1995)  for  modeling  full  elastic  wave  propagation  in  a  complex 
media. 

In  general,  the  boundary  integral  techniques  form  a  very  accurate  approach  to  modeling  the  propagation 
of  P  and  S  waves  in  an  irregularly  layered  media,  allowing  for  the  simple  inclusion  of  free  surface 
topography,  irregular  boundaries  separating  homogeneous  layers,  and  attenuation.  Unfortunately,  the 
boundary  integral  techniques  are  somewhat  limited  with  respect  to  regional  propagation,  since  the 
algorithm  quickly  becomes  computationally  intensive  as  the  boundary  conditions  along  each  interface 
must  be  solved  with  a  full  general  matrix  at  each  independent  frequency.  Bouchon  et  al.  (in  press)  have 
recently  shown  that  the  biconjugate  gradient  approach  combined  with  a  threshold  cutoff  in  the  matrix  can 
greatly  accelerate  these  algorithms.  However,  these  approaches  are  still  prohibitively  large  for  modeling 
high  frequency  responses  (/  >  1  Hz)  at  regional  distances.  We,  therefore,  implement  the  boundary 
integral  approach  more  as  a  site  tool  for  modeling  the  response  of  local  structure  to  incident  seismic 
energy.  In  addition,  this  code  has  proven  extremely  useful  for  benchmarking  the  validity  of  recently 
developed  wave  propagation  techniques  for  complex  media. 

Using  the  boundary  integral  technique  as  a  benchmarking  tool,  we  have  placed  a  major  emphasis  on 
making  finite-difference  an  accurate  and  viable  tool  for  modeling  wave  propagation  at  regional  distances. 
ELAS3D  has  become  a  highly  optimized  2D/3D  finite-difference  algorithm  which  is  fourth-order 
accurate  in  space  and  second  order  accurate  in  time.  This  code  has  a  run-time  visualization  feature  and 
post-processing  capabilities.  A  variable  density  grid  has  been  implemented,  yielding  significant  savings 
in  computer  speed  and  memory.  We  have  incorporated  both  active  grid  and  a  propagating  envelope 
algorithms.  These  speed  up  the  code  by  eliminating  calculations  in  regions  which  are  void  of 
“interesting”  seismic  energy.  In  addition,  we  have  implemented  a  very  efficient  scheme  for  modeling 
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the  response  of  free  surface  topography  in  the  2D  portion  of  the  code.  This  implementation  is  currently 
undergoing  rigorous  testing. 

We  have  ported  ELAS3D  to  several  platforms  including  the  Sun  workstations,  SGI,  IBM  2000,  and  the 
Cray  Y-MP,  giving  us  the  ability  to  model  wave  propagation  at  regional  distances.  Since  numerical 
dispersion  resulting  from  coarse  grid  samphng  is  the  greatest  limitation  to  propagation  at  larger 
distances,  we  have  started  implementing  ELAS3D  on  a  massively  parallel  processor  to  increase  the 
spatial  samphng  rate.  We  have  since  ported  a  simpler  acoustic  version  of  the  code  to  LLNL’s  256- 
processor  Meiko  CS-2  where  we  observed  a  speed  of  28  GFLOPS.  Initial  tests  show  that  the  parallelized 
version  of  the  full  elastic  code  will  dramatically  decrease  the  current  run-time  while  significantly 
increasing  the  available  memory.  The  largest  predicted  2D  model  is  300  x  3000  kilometers  at  10  Hz, 
given  a  constant  P-wave  velocity  of  5  km/s.  The  time  to  compute  a  full  seismogram  on  this  grid  is  on  the 
order  of  a  week.  Considering  a  typical  regional  propagation  distance  to  be  on  the  order  of  50  x  500  km, 
there  is  ample  memory  to  increase  the  sampling  rate  above  the  minimum  Emit  of  ten  points  per 
wavelength,  thereby  reducing  numerical  dispersion  considerably.  In  the  3D  case,  the  largest  possible 
model  is  25  x  25  x  25  kilometers  at  10  Hz  given  the  constant  P-wave  velocity  of  5  km/s.  The  time  to 
compute  a  whole  seismogram  is  on  the  order  of  only  a  few  hours. 

Application  to  Discriminant  variability 

An  improved  understanding  of  the  variability  of  regional  seismic  phases  with  distance  is  needed  to 
improve  the  performance  and  transportability  of  regional  seismic  discriminants.  Numerous  studies,  in  a 
number  of  regions,  have  observed  large  variations  in  each  of  the  dominant  regional  phases.  Studies  have 
shown  large  variations  in  Pn,  Sn,  and  Lg  over  relatively  short  distances  (e.g.  Keller  et  al.,  1994; 
Kadinsky-Cade  et  al.,  1981;  Ni  and  Barazangi,  1983).  An  improved  understanding  of  these  variations  has 
been  gained  from  numerous  empirical  studies  (e.g.  Chavez  and  Priestley,  1986;  Zhang  et  al.,  1994)  and 
theoretical  studies  (e.g.  Campillo,  1990;  Kennett,  1993).  The  numerical  approaches  discussed  in  the 
section  above,  give  us  the  tools  to  compliment  these  previous  studies  and  gain  physical  insight  into  the 
propagation  of  phases,  thereby  allowing  us  to  better  understand  the  transportability  of  discriminants. 

As  part  of  our  efforts  to  develop  procedures  for  regional  geophysical  characterization,  we  have  been 
making  and  assembling  empirical  observations  of  regional  phase  propagation.  At  the  same  time,  we  have 
been  interpreting  these  observations  with  our  numerical  modeling  to  describe  the  observed  phenomena. 
Examples  include  measurements  and  modehng  of  the  variability  of  regional  Pn,  Pg,  and  Lg  signals  from 
the  NPE  along  lines  to  the  east,  west,  and  northwest.  Figure  1  shows  the  three  arrays  deployed  to  record 
regional  phases  propagating  from  the  NPE.  The  330  km  long  Nevada  array  extends  from  the  NPE 
northward  through  the  Basin  and  Range  province  of  northwestern  Nevada.  The  3(X)  km  Arizona  array 
extends  eastward  from  the  Basin  and  Range  crossing  over  to  the  Colorado  plateau,  while  the  Sierra  array 
is  a  460  km  profile  which  extends  across  one  of  the  most  rugged  profiles  in  the  western  United  States. 
This  profile  crosses  Owen's  Valley,  the  Sierra  Nevada  Ranges,  the  Great  Valley,  and  finally  terminates 
just  short  of  the  San  Andreas  fault  zone.  The  average  station  separation  along  the  Sierra  line  is  less  than  2 
kilometers,  giving  excellent  station  coverage  of  regional  phases. 

Observations  and  modeling  of  variations  in  regional  phases,  both  with  distance  and  from  region  to 
region,  are  important  because  the  detection,  location,  and  identification  of  small  magnitude  events  rely 
heavily  on  observation  of  these  phases.  For  example,  one  of  the  most  popular  regional  discriminants  is 
based  on  the  ratio  of  high-frequency  (e.g.  6-8  Hz)  Pn  to  Lg  spectral  amplitudes.  If  one  of  these  phases  is 
more  sensitive  to  variations  in  structure  such  as  Moho  depth  or  topography,  we  would  expect  to  observe 
significant  variations  from  station  to  station,  or  at  a  single  station  for  events  from  different  azimuths. 
Figure  2  shows  recordings  of  the  NPE  along  the  line  that  extends  across  the  Sierra  Nevada  Range  and 
Figure  3a  demonstrates  the  variability  in  the  discriminants.  We  find  that  large  variations  in  Pn  map 
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directly  into  the  6-8  Hz  PnlLg  spectral  amplitude  ratio.  Variations  do  exist  in  Lg,  although  they  tend  to  be 
much  smaller  than  those  in  Pn.  If  such  variations  are  due  to  systematic  changes  in  features  like  crustal 
thickness  or  topography,  they  would  be  very  sensitive  to  source-receiver  distance  or  azimuth,  potentially 
introducing  large  variability  into  discriminants  such  as  PrdLg  spectral  ratios.  Based  on  the  SSCD 
refraction  experiment  (Ruppert  et  al.,  in  press;  Fliedner  et  al.,  in  press),  travel  time  tomography  has  been 
used  to  attain  a  3D  “ground  truth”  velocity  model  for  both  the  crust  and  mantle  in  this  region  using  data 
from  24  shots  and  three  profiles.  We  are  utilizing  this  profile  for  code  validation  as  we  try  to  better  define 
our  ability  to  predict  observed  variations  in  regional  phases.  This  testing  will  continue  and  improvements 
will  be  made  concurrent  with  our  efforts  of  regional  characterization  in  the  Middle  East  and  North  Africa. 

The  6-8  Hz  spectral  amplitude  of  the  three  most  dominant  regional  phases  is  plotted  as  a  function  of 
position  on  Figure  1.  There  is  a  substantial  amount  of  variation  in  the  amplitude  of  the  regional  phases, 
the  most  dramatic  variations  being  in  the  Pn  phase  while  both  the  Pg  and  Lg  phases  have  less  variation 
with  distance.  Although  variability  of  these  three  regional  phases  are  minimal  as  they  pass  through 
northwestern  Nevada,  there  is  an  order  of  magnitude  amplification  at  the  Broken  Hills  site  near  Fallon, 
Nevada.  Pn  is  amplified  by  almost  an  order  of  magnitude  relative  to  the  surrounding  stations,  while  Pg 
and  Lg  are  amplified  by  only  a  factor  of  two.  Figure  4  shows  two  proposed  structures  based  on  the  little 
that  we  know  of  the  geology  at  this  site  (Vitaliano,  1957)  and  the  observation  that  much  of  the 
amplification  consists  of  a  6  Hz  resonance.  Modeling  with  reflectivity  shows  that  although  the  spectral 
shape  and  the  drawn  out  time  domain  signal  can  be  generated  with  two  low  velocity  volcanic  layers,  one 
at  the  surface  and  the  other  a  few  hundred  meters  deep,  at  most  only  a  four  time  amplification  factor  is 
predicted.  Boundary  integral  modeling  of  a  basin  consisting  of  hydrothermally  altered  low-velocity 
volcanics  with  local  topography  demonstrates  that  2D  site  structure  can  extend  the  time  duration  of  the 
Pn  signal  and  provide  the  additional  amplification  observed  at  the  site. 

The  Arizona  line  also  shows  a  large  variation  in  Pn  with  distance,  while  Lg  again  tends  to  be  more  stable 
in  spectral  amplitude  at  6-8  Hz.  However,  the  oblique  polarization  of  Pn  at  the  site  of  amplification 
strongly  suggests  that  the  source  of  this  amplification  is  a  deeper  focusing  mechanism.  Zandt  (1995)  has 
shown  that  raytracing  combined  with  geologic  and  geophysical  control  along  this  profile  predict  crustal 
thickening  along  the  transition  between  the  Basin  and  Range  and  the  Colorado  Plateau,  resulting  in  two 
consecutive  Moho  steps  and  significant  lateral  velocity  variations.  This  model  predicts  focusing  at  the 
stations  where  Pn  appears  amplified  in  Figure  1.  Finite-difference  synthetics,  shown  in  Figure  5,  predict 
a  similar  amplification  of  Pn  at  the  corresponding  regional  distances.  Utilizing  the  broadband  nature  of 
this  finite-difference  modeling,  this  step  is  over  a  large  enough  distance  that,  in  addition  to  the  observed 
amplification  at  6  Hz,  significant  amplification  occurs  over  the  full  0.5  to  10  Hz  range  typically  utilized 
in  regional  discriminant  analysis. 

CONCLUSIONS  AND  RECOMMENDATIONS 

LLNL  will  continue  to  develop  a  numerical  foundation  for  the  synthesis  of  seismograms  at  regional 
distances.  We  have  developed  a  “ground  truth”  database  across  the  southern  Sierra  Nevadas  which  will 
be  used  for  the  continued  validation  of  our  numerical  modeling  capabilities.  In  addition,  future 
enhancements  will  be  made  to  ELAS3D  which  include  3D  implementations  of  attenuation,  topography, 
and  variable  density  grids,  along  with  a  2.5D  axial-symmetric  option.  This  code  will  be  made  available 
on  a  number  of  additional  platforms  including  workstation  clusters. 

Integrating  these  numerical  algorithms  with  our  empirical  study  of  phase  and  discriminant  variability  in 
the  southwestern  United  States,  we  have  demonstrated  the  great  variability  of  Pn  due  to  crustal 
thickening  and  local  site  resonances,  especieilly  as  observed  in  regions  containing  low-velocity  layers. 
Examples  include  the  crustal  thickening  of  the  Colorado  Plateau  and  Sierra  Nevada  range  and  near 
surface  low  velocity  structure  in  the  Great  Valley  and  Broken  Hills  region.  The  low  velocity  resonances 
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tend  to  amplify  Pn  both  due  to  the  resonance  itself  and  the  vertical  polarization  of  Pn.  Even  though  Lg 
is  amplified  by  these  low  velocities  it  becomes  radially  polarized,  significantly  reducing  the  amplitude  on 
the  vertical  component.  Therefore,  the  vertical  component  shows  an  anomalously  large  PnJLg  ratio, 
potentially  introducing  the  risk  of  false  alarms.  This  emphasizes  the  importance  of  three-component 
analysis  when  transporting  these  discriminants  to  different  regions.  One  example  where  numerical 
calculations  can  lead  to  insights  about  discriminant  performance  is  the  PnlLg  slope  discriminant 
(Goldstein,  1995).  This  discriminant  may  be  more  stable  than  taking  a  PnlLg  ratio  in  a  narrow  frequency 
band  because  of  the  inherent  averaging  that  occurs  over  narrow  band  resonances. 
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the  three  arrays  deployed  during  the  NPE.  The  three  pliase  amplitudes  are  plotted  side  by 
side  along  eaeh  array  so  that  Pn  is  plotted  over  the  aetual  station  location  and  Pg  and  Lg 
are  plotted  to  each  side  of  the  station  location.  Data  is  normalized  for  each  independent 
array  to  the  single  largest  phase  amplitude  recorded  at  a  station  along  the  array,  A  distance 
correction  has  not  been  applied  to  this  data. 
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Figure  2:  Plot  of  the  regional  seismograms  recorded  along  the  western  array.  Data  is  trace 
normalized  to  emphasize  phases  at  the  larger  offsets.  Lg  arrives  outside  the  time  window 
shown. 
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Figure  3;  Two  commonly  used  discriminants,  the  PnlLg  ratio  m  the  6-8  Hz  band  and  Pn! 
Lg  slope  over  the  full  spectral  range,  calculated  along  the  three  seismic  lines  deployed 
for  the  NPE.  We  have  applied  a  preliminary  distance  correction  to  this  data.  The  lower 
outliers  in  the  PnlLg  slope  are  partially  due  to  a  low  signal-to-noise  criterion. 
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Figure  5;  Finite-difference  synthetics  showing  the  response  of  the  Pn  and  Pg  phases  to 
crustal  thickening  along  the  transition  between  the  Basin  and  Range  and  the  Colorado 
Plateau.  At  200  km  the  Pn  is  amplified  by  approximately  a  factor  of  three  due  to  focusing 
from  a  Moho  step. 


